library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggpubr)
library(rstatix)
## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter
library(nparLD)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:rstatix':
## 
##     select
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(PMCMRplus)
library(WRS2)
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import numpy as np
from scipy.stats import pearsonr, spearmanr
import os
import gseapy as gp
from gseapy.parser import gsea_gmt_parser
from scipy.stats import wilcoxon, ttest_ind, ttest_rel
from statsmodels.stats.multitest import multipletests
from glob import glob
from statannotations.Annotator import Annotator
plt.rcParams.update({'font.size': 14})
annot = pd.read_csv('../data/batches12_chemo_regimen_annotation.tsv', sep='\t', index_col=0)
annot['Primary Study ID'] = annot.loc[:, 'Primary Study ID'].astype(str)
pt_reg = annot.set_index('Primary Study ID').loc[:, 'chemo_regimen'].to_dict()
platinum_type_dict = annot.set_index('Primary Study ID').loc[:, 'platinum_type'].to_dict()
def dyn_comp(u, ax, lin_dens3, log_dens3, title, xlabel, test='t-test_paired', test_name = 't-test paired', size=3, format_simple = True):
  sns.boxplot(y=u, x='sample_type', data = lin_dens3, ax=ax, order = ['TUR', 'Cystectomy'],
                    palette = tur_cys_palette, showfliers=False)
  sns.stripplot(y=u, x='sample_type', data = lin_dens3, ax=ax, order = ['TUR', 'Cystectomy'],
                    color='black', size=size)
      
  sns.lineplot(y=u, x='sample_type', data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.35, color = 'grey')
  ax.set_yscale("log")
  if format_simple:
    ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [('TUR', 'Cystectomy')], data = log_dens3, x='sample_type', y=u, order=['TUR', 'Cystectomy'])
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH')
  _, test_results = annotator.apply_test()._get_output()
  pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}\nN = {}'.format(test_name, pv, len(log_dens3.index))])
  annotator.annotate()
  print(u, 'Median TUR = {:.3g}'.format(lin_dens3[lin_dens3.sample_type == 'TUR'].loc[:, u].median()),
  '\nMedian Cystectomy = {:.3g}'.format(lin_dens3[lin_dens3.sample_type == 'Cystectomy'].loc[:, u].median()))
  
def dyn_comp_lin(u, ax, lin_dens3, title, xlabel, test='t-test_paired', test_name = 't-test paired', size=3):
  sns.boxplot(y=u, x='sample_type', data = lin_dens3, ax=ax, order = ['TUR', 'Cystectomy'],
                    palette = tur_cys_palette, showfliers=False)
  sns.stripplot(y=u, x='sample_type', data = lin_dens3, ax=ax, order = ['TUR', 'Cystectomy'],
                    color='black', size=size)
      
  sns.lineplot(y=u, x='sample_type', data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.35, color = 'grey')
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [('TUR', 'Cystectomy')], data = lin_dens3, x='sample_type', y=u, order=['TUR', 'Cystectomy'])
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH')
  _, test_results = annotator.apply_test()._get_output()
  pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}\nN = {}'.format(test_name, pv, len(lin_dens3.index))])
  annotator.annotate()
  print(u, 'Median TUR = {:.3g}'.format(lin_dens3[lin_dens3.sample_type == 'TUR'].loc[:, u].median()),
  '\nMedian Cystectomy = {:.3g}'.format(lin_dens3[lin_dens3.sample_type == 'Cystectomy'].loc[:, u].median()))
  
def dyn_comp_delta(u, ax, lin_dens3, title, xlabel, to_compare, order, test = 't-test_paired', size=3, ofs = 1):
  sns.boxplot(y=u, x=to_compare, data = lin_dens3, ax=ax, order = order,
                    color = 'mediumseagreen', showfliers=False)
  sns.swarmplot(y=u, x=to_compare, data = lin_dens3, ax=ax, order = order,
                    color='black', size=size)
      
  sns.lineplot(y=u, x=to_compare, data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.7)
      
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [order], data = lin_dens3, x=to_compare, y=u, order=order)
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH', line_offset = ofs, line_offset_to_group = ofs)
  _, test_results = annotator.apply_test()._get_output()
  pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}'.format(test, pv)])
  annotator.annotate(line_offset = ofs, line_offset_to_group = ofs)

Figure 2

norm_log_dens = pd.read_csv('../data/combined_norm_log10_density.tsv', sep='\t', index_col=0)
norm_log_dens_tumor = pd.read_csv('../data/tumor_norm_log10_density.tsv', sep='\t', index_col=0)
norm_log_dens_stroma = pd.read_csv('../data/stroma_norm_log10_density.tsv', sep='\t', index_col=0)
raw_dens = pd.read_csv('../data/all_densities_combined_kde_mm2.tsv', sep='\t', index_col=0)
dens_to_comp = ['B_cells', 'CD4_Tcells', 'Macrophages', 'Tregs', 'CD8_Tcells', 'PanCK+_cells', 'Negative_cells']
def transform_dens(norm_log_dens, addition = 0):
  norm_log_dens = norm_log_dens.loc[:, dens_to_comp] + addition
  lin_dens2 = (10**norm_log_dens).assign(sample_type = norm_log_dens.index.str.split(' ').str[1]).assign(pt_number = norm_log_dens.index.str.split(' ').str[0])
  lin_dens3 = lin_dens2[lin_dens2.pt_number.isin(lin_dens2.pt_number.value_counts()[lin_dens2.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False)
  lin_dens3 = lin_dens3.assign(chemo_regimen = lin_dens3.pt_number.map(pt_reg)).assign(platinum_type = lin_dens3.pt_number.map(platinum_type_dict))
  log_dens2 = norm_log_dens.assign(sample_type = norm_log_dens.index.str.split(' ').str[1]).assign(pt_number = norm_log_dens.index.str.split(' ').str[0])
  log_dens3 = log_dens2[log_dens2.pt_number.isin(log_dens2.pt_number.value_counts()[log_dens2.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False)
  log_dens3 = log_dens3.assign(chemo_regimen = log_dens3.pt_number.map(pt_reg)).assign(platinum_type = log_dens3.pt_number.map(platinum_type_dict))
  return lin_dens3, log_dens3
lin_dens3, log_dens3 = transform_dens(norm_log_dens, addition = 2)
raw_lin_dens3, raw_log_dens3 = transform_dens(raw_dens)
lin_dens3_tumor, log_dens3_tumor = transform_dens(norm_log_dens_tumor, addition = 2)
lin_dens3_stroma, log_dens3_stroma = transform_dens(norm_log_dens_stroma, addition = 2)
pdl1 = pd.read_csv('../data/Prepost NAC patients overview PD-L1.csv', sep=';', index_col=0)
pdl1_cys = pdl1.loc[~pdl1.index.isna(), ['Material.1', 'CPS (%).1', 'IC (%).1', 'TC (%).1']].dropna()
pdl1_tur = pdl1.loc[~pdl1.index.isna(), ['Material', 'CPS (%)', 'IC (%)', 'TC (%)']].dropna().replace('76,5', '76.5')
pdl1_tur.index = pdl1_tur.index.astype(int).astype(str) + ' TUR'
pdl1_cys.index = pdl1_cys.index.astype(int).astype(str) + ' Cystectomy'
pdl1_cys.columns = ['Material', 'CPS (%)', 'IC (%)', 'TC (%)']
pdl1_tur = pdl1_tur.drop('Material', axis=1).astype(float)+1
pdl1_cys = pdl1_cys.drop('Material', axis=1).astype(float)+1
pdl1_df = pd.concat([pdl1_tur.assign(sample_type = 'TUR'), pdl1_cys.assign(sample_type = 'Cystectomy')])
pdl1_df_log = pd.concat([np.log10(pdl1_tur).assign(sample_type = 'TUR'), np.log10(pdl1_cys).assign(sample_type = 'Cystectomy')])
pdl1_df2 = pdl1_df.assign(pt_number = pdl1_df.index.str.split(' ').str[0])
pdl1_df3 = pdl1_df2[pdl1_df2.pt_number.isin(pdl1_df2.pt_number.value_counts()[pdl1_df2.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False)
pdl1_df3 = pdl1_df3.assign(platinum_type = pdl1_df3.pt_number.map(platinum_type_dict)).assign(chemo_regimen = pdl1_df3.pt_number.map(pt_reg))

pdl1_df2_log = pdl1_df_log.assign(pt_number = pdl1_df_log.index.str.split(' ').str[0])
pdl1_df3_log = pdl1_df2_log[pdl1_df2_log.pt_number.isin(pdl1_df2_log.pt_number.value_counts()[pdl1_df2_log.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False)
pdl1_df3_log = pdl1_df3_log.assign(platinum_type = pdl1_df3_log.pt_number.map(platinum_type_dict)).assign(chemo_regimen = pdl1_df3_log.pt_number.map(pt_reg))
ssgsea = pd.read_csv('../data/ssgsea_scores_Fig2.tsv', sep='\t', index_col=0)
ssgsea = ssgsea.reset_index().assign(sample_type = ssgsea.reset_index().loc[:, 'index'].apply(lambda x: x.split(' ')[1])).assign(pt_number = ssgsea.reset_index().loc[:, 'index'].apply(lambda x: x.split(' ')[0])).set_index('index')
ps_ssgsea0 = ssgsea[ssgsea.pt_number.isin(ssgsea.pt_number.value_counts()[ssgsea.pt_number.value_counts() == 2].index)].sort_values(by = ['pt_number', 'sample_type'], ascending=False)
ps_ssgsea0 = ps_ssgsea0.assign(chemo_regimen = ps_ssgsea0.pt_number.map(pt_reg)).assign(platinum_type = ps_ssgsea0.pt_number.map(platinum_type_dict))
tur_cys_palette = {'TUR': 'lightcoral', 'Cystectomy': 'lightskyblue'}
plt.close()
fig, axs = plt.subplots(3, 3, figsize = (15, 15))
af = axs.flat
dyn_comp('CD8_Tcells', af[0], lin_dens3_tumor,  log_dens3_tumor, '', 'CD8 T cell percentage in tumor, %')
## TUR vs. Cystectomy: t-test paired p=0.015
## N = 170
## CD8_Tcells Median TUR = 1.33 
## Median Cystectomy = 1.95
dyn_comp('CD8_Tcells', af[1], lin_dens3_stroma,  log_dens3_stroma, '', 'CD8 T cell percentage in stroma, %')
## TUR vs. Cystectomy: t-test paired p=0.017
## N = 170
## CD8_Tcells Median TUR = 5.1 
## Median Cystectomy = 6.51
dyn_comp('CD8_Tcells', af[2], lin_dens3,  log_dens3, '', 'CD8 T cell percentage in tumor + stroma, %')
## TUR vs. Cystectomy: t-test paired p=0.00076
## N = 170
## CD8_Tcells Median TUR = 2.34 
## Median Cystectomy = 4.17
_ = af[0].set_ylim((0.01), 400)
_ = af[1].set_ylim((0.01), 400)
_ = af[2].set_ylim((0.01), 400)
for i in range(3):
  _ = af[i].set_yticks([0.01,0.1,1, 10, 100], ['0.01','0.1','1', '10', '100'])

for i,u in enumerate(['IC (%)', 'TC (%)', 'CPS (%)']):
  dyn_comp(u, af[i+3], pdl1_df3, pdl1_df3_log, '', '{} PDL1 score'.format(u), test='Wilcoxon', test_name='Wilcoxon')
  _ = af[i+3].set_ylim((0.5, 400))
  _ = af[i+3].set_yticks([1,10,100], ['1','10', '100'])
  

ps_ssgsea = ps_ssgsea0.copy()
dyn_comp_lin('TGFB_Mariathasan', af[6], ps_ssgsea, '', 'fibroblast TGFb ssGSEA score')
dyn_comp_lin('IFNG_signature', af[7], ps_ssgsea, '', 'INFg ssGSEA score')
dyn_comp_lin('CD8_T_EFFECTOR', af[8], ps_ssgsea, '', 'effector CD8 T cell ssGSEA score')
for i in range(9):
  _ = af[i].set_xticks(['TUR', 'Cystectomy'], ['TUR-BT', 'Cystectomy'])
plt.subplots_adjust(hspace=0.2, wspace=0.4)
plt.show()

# Figure S3

plt.close()
fig, axs = plt.subplots(2, 6, figsize = (22.7, 11))
af = axs.flat
cell_names = {'CD8_Tcells': 'CD8 T cell', 'Macrophages': 'Macrophage', 'CD4_Tcells': 'T helper',  'Tregs': 'Treg', 'B_cells': 'B cell', 'PanCK+_cells': 'Cancer cell', 'Negative_cells': 'Negative cell'}
for i, u in enumerate(['CD8_Tcells', 'Macrophages', 'CD4_Tcells',  'Tregs', 'B_cells', 'PanCK+_cells']):
  dyn_comp(u, af[i], lin_dens3,  log_dens3, '', '{} percentage, %'.format(cell_names[u]), test='Wilcoxon', test_name='Wilcoxon')
  _ = af[i].set_xticks(['TUR', 'Cystectomy'], ['TUR-BT', 'Cystectomy'])
  
for i, u in enumerate(['CD8_Tcells', 'Macrophages', 'CD4_Tcells',  'Tregs', 'B_cells', 'PanCK+_cells']):
  dyn_comp(u, af[i+6], raw_lin_dens3,  raw_log_dens3, '', '{} density, cells/mm2'.format(cell_names[u]), test='Wilcoxon', test_name='Wilcoxon', format_simple = False)
  _ = af[i+6].set_xticks(['TUR', 'Cystectomy'], ['TUR-BT', 'Cystectomy'])

plt.subplots_adjust(hspace=0.2, wspace=0.47)
plt.show()

cys_den = log_dens3[log_dens3.sample_type == 'Cystectomy'].set_index('pt_number')
tur_den = log_dens3[log_dens3.sample_type == 'TUR'].set_index('pt_number')
common_in = cys_den.index.intersection(tur_den.index)
delta_den = cys_den.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1) - tur_den.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1)
delta_den = delta_den.assign(pt_number = delta_den.index)
delta_den = delta_den.assign(chemo_regimen = delta_den.pt_number.map(pt_reg)).assign(platinum_type = delta_den.pt_number.map(platinum_type_dict))
cys_ps= ps_ssgsea0[ps_ssgsea0.sample_type == 'Cystectomy'].set_index('pt_number')
tur_ps= ps_ssgsea0[ps_ssgsea0.sample_type == 'TUR'].set_index('pt_number')
common_in = cys_ps.index.intersection(tur_ps.index)
delta_ps = cys_ps.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1) - tur_ps.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1)
delta_ps = delta_ps.assign(pt_number = delta_ps.index)
delta_ps = delta_ps.assign(chemo_regimen = delta_ps.pt_number.map(pt_reg)).assign(platinum_type = delta_ps.pt_number.map(platinum_type_dict))

Figure 3D

all_med = pd.read_csv('../data/all_median_distances_wo_offset.tsv', sep='\t', index_col=0)
all_med_log = np.log10(all_med.drop(['sample_type', 'pt_number'], axis=1)).assign(sample_type = all_med.index.str.split(' ').str[1]).assign(pt_number = all_med.index.str.split(' ').str[0])
all_med_lin = all_med.drop(['sample_type'], axis=1).assign(sample_type = all_med.index.str.split(' ').str[1]).assign(pt_number = all_med.index.str.split(' ').str[0])
all_med_log = all_med_log[all_med_log.pt_number.isin(all_med_log.pt_number.value_counts()[all_med_log.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False).assign(chemo_regimen = all_med_log.pt_number.map(pt_reg)).assign(platinum_type =all_med_log.pt_number.map(platinum_type_dict))
all_med_lin = all_med_lin[all_med_lin.pt_number.isin(all_med_lin.pt_number.value_counts()[all_med_lin.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False).assign(chemo_regimen = all_med_lin.pt_number.map(pt_reg)).assign(platinum_type = all_med_lin.pt_number.map(platinum_type_dict))
plt.close()
fig, axs = plt.subplots(1, 2, figsize = (8, 5))
af = axs.flat
dyn_comp('CD8_Tcell_to_tumor_median', af[0], all_med_lin,  all_med_log, '', 'CD8 T cell to 1-NN cancer cell\nmedian distance, um', test='Wilcoxon', test_name='Wilcoxon')
dyn_comp('Macrophage_to_tumor_median', af[1], all_med_lin,  all_med_log, '', 'Macrophage to 1-NN cancer cell\nmedian distance, um', test='Wilcoxon', test_name='Wilcoxon')
af[0].set_yticks([10, 50, 100, 500], ['10','50','100','500'])
af[1].set_yticks([10, 50, 100], ['10','50','100'])
#af[0].set_ylim((0.0001), 2)
#af[1].set_ylim((0.0001), 2)
for i in range(2):
  af[i].set_xticks(['TUR', 'Cystectomy'], ['TUR-BT', 'Cystectomy'])
#_ = fig.set_tight_layout(True)
plt.subplots_adjust(hspace=0.2, wspace=0.5)
plt.show()

Figure 5

tur_dist2 = (all_med_log[all_med_log.sample_type == 'TUR']).set_index('pt_number')
cys_dist2 = (all_med_log[all_med_log.sample_type == 'Cystectomy']).set_index('pt_number')
common_index0 = tur_dist2.index.intersection(cys_dist2.index)
delta_dist0 = cys_dist2.loc[common_index0, ['CD8_Tcell_to_tumor_median', 'Macrophage_to_tumor_median']] - tur_dist2.loc[common_index0, ['CD8_Tcell_to_tumor_median', 'Macrophage_to_tumor_median']]
delta_dist0 = delta_dist0.assign(pt_number = delta_dist0.index.str.split(' ').str[0])
delta_dist0 = delta_dist0.assign(chemo_regimen = delta_dist0.pt_number.map(pt_reg)).assign(platinum_type = delta_dist0.pt_number.map(platinum_type_dict)).dropna(subset = ['chemo_regimen'])
pdl1_df3_log = pdl1_df3_log.rename(columns = {'IC (%)': 'IC', 'CPS (%)': 'CPS', 'TC (%)': 'TC'}).dropna()
pdl1_df3 = pdl1_df3.rename(columns = {'IC (%)': 'IC', 'CPS (%)': 'CPS', 'TC (%)': 'TC'}).dropna()
cys_pdl = pdl1_df3_log[pdl1_df3_log.sample_type == 'Cystectomy'].set_index('pt_number')
tur_pdl = pdl1_df3_log[pdl1_df3_log.sample_type == 'TUR'].set_index('pt_number')
common_in = cys_pdl.index.intersection(tur_pdl.index)
delta_pdl = cys_pdl.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1) - tur_pdl.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1)
delta_pdl = delta_pdl.assign(pt_number = delta_pdl.index)
delta_pdl = delta_pdl.assign(chemo_regimen = delta_pdl.pt_number.map(pt_reg)).assign(platinum_type = delta_pdl.pt_number.map(platinum_type_dict))
def dyn_comp_anova(u, ax, lin_dens3, log_dens3, title, xlabel, loc_leg, test='t-test_paired', test_name='t-test paired'):
  sns.boxplot(y=u, x = 'chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'],
                    palette = tur_cys_palette, showfliers=False)
  sns.stripplot(y=u, x = 'chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'],
                    color='black', dodge=True, size=3)
      
  #sns.lineplot(y=u, x='chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.5, color = 'grey')
  ax.set_yscale("log")    
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [(('MVAC', 'TUR'), ('MVAC', 'Cystectomy')), (('gemcitabine_platinum', 'TUR'), ('gemcitabine_platinum', 'Cystectomy'))], data = log_dens3, x='chemo_regimen', hue = 'sample_type', y=u, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'])
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH')
  _, test_results = annotator.apply_test()._get_output()
  #pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}'.format(test_name, test_results[0].data.pvalue), '{} p={:.2g}'.format(test_name, test_results[1].data.pvalue)])
  annotator.annotate()
  handles, labels = ax.get_legend_handles_labels()
  ax.legend(handles[0:2], ['TUR-BT', 'Cystectomy'], loc=loc_leg)
  
def dyn_comp_lin_anova(u, ax, lin_dens3, title, xlabel, loc_leg, test='t-test_paired', test_name='t-test paired'):
  sns.boxplot(y=u, x = 'chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'],
                    palette = tur_cys_palette, showfliers=False)
  sns.stripplot(y=u, x = 'chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'],
                    color='black', dodge=True, size=3)
      
  #sns.lineplot(y=u, hue='chemo_regimen', x='sample_type', data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.5, color = 'grey')
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [(('MVAC', 'TUR'), ('MVAC', 'Cystectomy')), (('gemcitabine_platinum', 'TUR'), ('gemcitabine_platinum', 'Cystectomy'))], data = lin_dens3, x='chemo_regimen', hue = 'sample_type', y=u, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'])
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH')
  _, test_results = annotator.apply_test()._get_output()
  #pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}'.format(test_name, test_results[0].data.pvalue), '{} paired p={:.2g}'.format(test_name, test_results[1].data.pvalue)])
  annotator.annotate()
  handles, labels = ax.get_legend_handles_labels()
  ax.legend(handles[0:2], ['TUR-BT', 'Cystectomy'], loc=loc_leg)
plt.rcParams.update({'font.size': 15})
testt = 'Wilcoxon'
testt_name = 'Wilcoxon'
plt.close()
fig, axs = plt.subplots(2, 2, figsize = (14.5, 12.5), width_ratios = [10, 5])
af = axs.flat
u = 'Macrophages'
dyn_comp_anova(u, af[0], lin_dens3, log_dens3, '', 'Macrophage cell percentage', loc_leg=(0.3, 0.03), test='Wilcoxon', test_name = 'Wilcoxon')
dyn_comp_delta(u, af[1], delta_den, '', 'Δ log10(macrophage percentage)\nupon treatment', to_compare = 'chemo_regimen', order = ('MVAC', 'gemcitabine_platinum'), test = 'Mann-Whitney', ofs = -0.08, size=5)
af[0].set_xticklabels(['MVAC', 'gemcitabine +\nplatinum'])
af[1].set_xticklabels(['MVAC', 'gemcitabine +\nplatinum'])
af[1].set_xlim((-0.7,1.7))
af[1].set_ylim((-1.4,2.8))
af[0].set_yticks([0.1,1, 10, 100], ['0.1','1', '10', '100'])
u = 'TGFB_Mariathasan'
dyn_comp_lin_anova(u, af[2], ps_ssgsea0, '', 'Fibroblast TGFb ssGSEA score', loc_leg='lower right', test='Wilcoxon', test_name = 'Wilcoxon')
dyn_comp_delta(u, af[3], delta_ps, '', 'Δ fibroblast TGFb upon treatment', to_compare = 'chemo_regimen', order = ('MVAC', 'gemcitabine_platinum'), test = 'Mann-Whitney', ofs=0.1, size=5)
af[2].set_xticklabels(['MVAC', 'gemcitabine +\nplatinum'])
af[3].set_xticklabels(['MVAC', 'gemcitabine +\nplatinum'])
af[3].set_xlim((-0.7,1.7))
#_ = fig.set_tight_layout(True)
plt.subplots_adjust(hspace=0.3, wspace=0.25)
plt.show()

Mixed ANOVA models

densities

lin_dens3_cr = lin_dens3[lin_dens3.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]
#lin_dens3_ptt = lin_dens3[lin_dens3.platinum_type.isin(['Cisplatin', 'Carboplatin'])]
log_dens3_cr = log_dens3[log_dens3.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]
log_dens3_ptt = log_dens3[log_dens3.platinum_type.isin(['Cisplatin', 'Carboplatin'])]

CD8 T cells

py$log_dens3_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(CD8_Tcells)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable   statistic     p
##   <chr>       <chr>                <chr>          <dbl> <dbl>
## 1 Cystectomy  MVAC                 CD8_Tcells     0.985 0.936
## 2 Cystectomy  gemcitabine_platinum CD8_Tcells     0.984 0.745
## 3 TUR         MVAC                 CD8_Tcells     0.969 0.499
## 4 TUR         gemcitabine_platinum CD8_Tcells     0.984 0.759
py$log_dens3_ptt %>%
  group_by(sample_type, platinum_type) %>%
  shapiro_test(CD8_Tcells)
## # A tibble: 4 × 5
##   sample_type platinum_type variable   statistic     p
##   <chr>       <chr>         <chr>          <dbl> <dbl>
## 1 Cystectomy  Carboplatin   CD8_Tcells     0.910 0.351
## 2 Cystectomy  Cisplatin     CD8_Tcells     0.991 0.909
## 3 TUR         Carboplatin   CD8_Tcells     0.923 0.454
## 4 TUR         Cisplatin     CD8_Tcells     0.981 0.370
py$log_dens3_cr %>%
  group_by(sample_type) %>%
  levene_test(CD8_Tcells ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    78     0.296 0.588
## 2 TUR             1    78     1.26  0.264
py$log_dens3_ptt %>%
  group_by(sample_type) %>%
  levene_test(CD8_Tcells ~ platinum_type)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic      p
##   <chr>       <int> <int>     <dbl>  <dbl>
## 1 Cystectomy      1    74     0.753 0.388 
## 2 TUR             1    74     2.88  0.0939
box_m(py$log_dens3_cr[, "CD8_Tcells", drop = FALSE], py$log_dens3_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      1.10   0.294         1 Box's M-test for Homogeneity of Covariance Matric…
box_m(py$log_dens3_ptt[, "CD8_Tcells", drop = FALSE], py$log_dens3_ptt$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      1.22   0.270         1 Box's M-test for Homogeneity of Covariance Matric…
res.aov <- anova_test(
  data = py$log_dens3_cr, dv = CD8_Tcells, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F        p p<.05      ges
## 1             chemo_regimen   1  78  0.048 0.827000       0.000431
## 2               sample_type   1  78 13.611 0.000414     * 0.050000
## 3 chemo_regimen:sample_type   1  78  0.857 0.358000       0.003000
bwtrim(CD8_Tcells ~ chemo_regimen*sample_type, id = pt_number, data = py$log_dens3_cr)
## Call:
## bwtrim(formula = CD8_Tcells ~ chemo_regimen * sample_type, id = pt_number, 
##     data = py$log_dens3_cr)
## 
##                            value df1     df2 p.value
## chemo_regimen             0.0381   1 56.2559  0.8459
## sample_type               8.3784   1 49.8259  0.0056
## chemo_regimen:sample_type 0.2000   1 49.8259  0.6566
res.aov <- anova_test(
  data = py$log_dens3_ptt, dv = CD8_Tcells, wid = pt_number,
  between = platinum_type, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F     p p<.05   ges
## 1             platinum_type   1  74  1.591 0.211       0.015
## 2               sample_type   1  74 10.900 0.001     * 0.043
## 3 platinum_type:sample_type   1  74  1.302 0.258       0.005
bwtrim(CD8_Tcells ~ platinum_type*sample_type, id = pt_number, data = py$log_dens3_ptt)
## Call:
## bwtrim(formula = CD8_Tcells ~ platinum_type * sample_type, id = pt_number, 
##     data = py$log_dens3_ptt)
## 
##                             value df1    df2 p.value
## platinum_type              3.8513   1 7.9619  0.0855
## sample_type               13.2831   1 6.5957  0.0091
## platinum_type:sample_type  3.1792   1 6.5957  0.1204

Macrophages

py$log_dens3_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(Macrophages)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable    statistic         p
##   <chr>       <chr>                <chr>           <dbl>     <dbl>
## 1 Cystectomy  MVAC                 Macrophages     0.969 0.499    
## 2 Cystectomy  gemcitabine_platinum Macrophages     0.876 0.0000984
## 3 TUR         MVAC                 Macrophages     0.881 0.00254  
## 4 TUR         gemcitabine_platinum Macrophages     0.957 0.0695
py$log_dens3_ptt %>%
  group_by(sample_type, platinum_type) %>%
  shapiro_test(Macrophages)
## # A tibble: 4 × 5
##   sample_type platinum_type variable    statistic         p
##   <chr>       <chr>         <chr>           <dbl>     <dbl>
## 1 Cystectomy  Carboplatin   Macrophages     0.923 0.458    
## 2 Cystectomy  Cisplatin     Macrophages     0.902 0.0000603
## 3 TUR         Carboplatin   Macrophages     0.933 0.541    
## 4 TUR         Cisplatin     Macrophages     0.926 0.000588
py$log_dens3_cr %>%
  group_by(sample_type) %>%
  levene_test(Macrophages ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    78    0.929  0.338
## 2 TUR             1    78    0.0228 0.880
box_m(py$log_dens3_cr[, "Macrophages", drop = FALSE], py$log_dens3_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1     0.812   0.367         1 Box's M-test for Homogeneity of Covariance Matric…
bwtrim(Macrophages ~ chemo_regimen*sample_type, id = pt_number, data = py$log_dens3_cr)
## Call:
## bwtrim(formula = Macrophages ~ chemo_regimen * sample_type, id = pt_number, 
##     data = py$log_dens3_cr)
## 
##                             value df1     df2 p.value
## chemo_regimen              2.8689   1 59.9896  0.0955
## sample_type               13.0202   1 55.4999  0.0007
## chemo_regimen:sample_type  2.2379   1 55.4999  0.1403
res.aov <- anova_test(
  data = py$log_dens3_cr, dv = Macrophages, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F        p p<.05   ges
## 1             chemo_regimen   1  78  2.467 0.120000       0.021
## 2               sample_type   1  78 14.688 0.000255     * 0.056
## 3 chemo_regimen:sample_type   1  78  5.646 0.020000     * 0.022
f1.ld.f1(y = py$log_dens3_cr$Macrophages, time = py$log_dens3_cr$sample_type, group = py$log_dens3_cr$chemo_regimen, subject = py$log_dens3_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##            Statistic df      p-value
## Group       3.184968  1 0.0743184625
## Time       14.764237  1 0.0001218242
## Group:Time  5.140979  1 0.0233674643
f1.ld.f1(y = py$log_dens3_ptt$Macrophages, time = py$log_dens3_ptt$sample_type, group = py$log_dens3_ptt$platinum_type, subject = py$log_dens3_ptt$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('Cisplatin', 'Carboplatin'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   Cisplatin Carboplatin 
##  If the order is not correct, specify the correct order in time.order or group.order.

##            Statistic df     p-value
## Group      9.6853793  1 0.001857400
## Time       7.2927779  1 0.006923234
## Group:Time 0.2585108  1 0.611145107

TGFb signature

ps_ssgsea_cr = ps_ssgsea0[ps_ssgsea0.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]#.drop(['15 TUR', '15 Cystectomy', '7 TUR', '7 Cystectomy', '31 TUR', '31 Cystectomy'])
ps_ssgsea_ptt = ps_ssgsea0[ps_ssgsea0.platinum_type.isin(['Cisplatin', 'Carboplatin'])]
py$ps_ssgsea_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(TGFB_Mariathasan)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable         statistic     p
##   <chr>       <chr>                <chr>                <dbl> <dbl>
## 1 Cystectomy  MVAC                 TGFB_Mariathasan     0.936 0.164
## 2 Cystectomy  gemcitabine_platinum TGFB_Mariathasan     0.964 0.144
## 3 TUR         MVAC                 TGFB_Mariathasan     0.926 0.102
## 4 TUR         gemcitabine_platinum TGFB_Mariathasan     0.989 0.921
py$ps_ssgsea_ptt %>%
  group_by(sample_type, platinum_type) %>%
  shapiro_test(TGFB_Mariathasan)
## # A tibble: 4 × 5
##   sample_type platinum_type variable         statistic      p
##   <chr>       <chr>         <chr>                <dbl>  <dbl>
## 1 Cystectomy  Carboplatin   TGFB_Mariathasan     0.889 0.134 
## 2 Cystectomy  Cisplatin     TGFB_Mariathasan     0.957 0.0427
## 3 TUR         Carboplatin   TGFB_Mariathasan     0.896 0.167 
## 4 TUR         Cisplatin     TGFB_Mariathasan     0.977 0.361
py$ps_ssgsea_cr %>%
  group_by(sample_type) %>%
  levene_test(TGFB_Mariathasan ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic      p
##   <chr>       <int> <int>     <dbl>  <dbl>
## 1 Cystectomy      1    68     0.454 0.503 
## 2 TUR             1    68     6.11  0.0159
py$ps_ssgsea_ptt %>%
  group_by(sample_type) %>%
  levene_test(TGFB_Mariathasan ~ platinum_type)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    66    0.777  0.381
## 2 TUR             1    66    0.0422 0.838
box_m(py$ps_ssgsea_cr[, "TGFB_Mariathasan", drop = FALSE], py$ps_ssgsea_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      3.47  0.0627         1 Box's M-test for Homogeneity of Covariance Matric…
box_m(py$ps_ssgsea_ptt[, "TGFB_Mariathasan", drop = FALSE], py$ps_ssgsea_ptt$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      3.36  0.0668         1 Box's M-test for Homogeneity of Covariance Matric…
res.aov <- anova_test(
  data = py$ps_ssgsea_cr, dv = TGFB_Mariathasan, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F        p p<.05   ges
## 1             chemo_regimen   1  68  1.101 2.98e-01       0.010
## 2               sample_type   1  68 24.729 4.74e-06     * 0.123
## 3 chemo_regimen:sample_type   1  68  4.629 3.50e-02     * 0.026
tgf_df = bwtrim(TGFB_Mariathasan ~ chemo_regimen*sample_type, id = pt_number, data = py$ps_ssgsea_cr)
bwtrim(TGFB_Mariathasan ~ platinum_type*sample_type, id = pt_number, data = py$ps_ssgsea_ptt)
## Call:
## bwtrim(formula = TGFB_Mariathasan ~ platinum_type * sample_type, 
##     id = pt_number, data = py$ps_ssgsea_ptt)
## 
##                             value df1    df2 p.value
## platinum_type              0.1571   1 8.7079  0.7014
## sample_type               26.5357   1 9.1856  0.0006
## platinum_type:sample_type  0.9076   1 9.1856  0.3652
f1.ld.f1(y = py$ps_ssgsea_cr$TGFB_Mariathasan, time = py$ps_ssgsea_cr$sample_type, group = py$ps_ssgsea_cr$chemo_regimen, subject = py$ps_ssgsea_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##             Statistic df      p-value
## Group       0.4801879  1 4.883372e-01
## Time       27.8105919  1 1.337916e-07
## Group:Time  3.1356103  1 7.659971e-02
f1.ld.f1(y = py$ps_ssgsea_ptt$TGFB_Mariathasan, time = py$ps_ssgsea_ptt$sample_type, group = py$ps_ssgsea_ptt$platinum_type, subject = py$ps_ssgsea_ptt$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('Cisplatin', 'Carboplatin'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   Cisplatin Carboplatin 
##  If the order is not correct, specify the correct order in time.order or group.order.

##             Statistic df      p-value
## Group       0.1589382  1 6.901360e-01
## Time       34.3513238  1 4.600873e-09
## Group:Time  1.7966710  1 1.801155e-01

PDL1 score

pdl_cr = pdl1_df3_log[pdl1_df3_log.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]
pdl_ptt = pdl1_df3_log[pdl1_df3_log.platinum_type.isin(['Cisplatin', 'Carboplatin'])]
py$pdl_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(IC)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable statistic         p
##   <chr>       <chr>                <chr>        <dbl>     <dbl>
## 1 Cystectomy  MVAC                 IC           0.797 0.0000351
## 2 Cystectomy  gemcitabine_platinum IC           0.883 0.000454 
## 3 TUR         MVAC                 IC           0.867 0.000994 
## 4 TUR         gemcitabine_platinum IC           0.904 0.00184
py$pdl_ptt %>%
  group_by(sample_type, platinum_type) %>%
  shapiro_test(IC)
## # A tibble: 4 × 5
##   sample_type platinum_type variable statistic           p
##   <chr>       <chr>         <chr>        <dbl>       <dbl>
## 1 Cystectomy  Carboplatin   IC           0.932 0.570      
## 2 Cystectomy  Cisplatin     IC           0.814 0.000000182
## 3 TUR         Carboplatin   IC           0.824 0.0698     
## 4 TUR         Cisplatin     IC           0.895 0.0000616
py$pdl_cr %>%
  group_by(sample_type) %>%
  levene_test(IC ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    72    2.36   0.129
## 2 TUR             1    72    0.0119 0.913
box_m(py$pdl_cr[, "IC", drop = FALSE], py$pdl_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      1.34   0.246         1 Box's M-test for Homogeneity of Covariance Matric…
#bwtrim(IC ~ chemo_regimen*sample_type, id = pt_number, data = py$pdl_cr)
#bwtrim(IC ~ chemo_regimen*sample_type, id = pt_number, data = py$pdl_ptt)
f1.ld.f1(y = py$pdl_cr$IC, time = py$pdl_cr$sample_type, group = py$pdl_cr$chemo_regimen, subject = py$pdl_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##            Statistic df      p-value
## Group       6.125007  1 1.332828e-02
## Time       19.517389  1 9.968807e-06
## Group:Time  2.295293  1 1.297667e-01
f1.ld.f1(y = py$pdl_ptt$IC, time = py$pdl_ptt$sample_type, group = py$pdl_ptt$platinum_type, subject = py$pdl_ptt$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('Cisplatin', 'Carboplatin'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   Cisplatin Carboplatin 
##  If the order is not correct, specify the correct order in time.order or group.order.

##              Statistic df      p-value
## Group       0.06227595  1 8.029342e-01
## Time       22.36575350  1 2.253569e-06
## Group:Time  0.72560845  1 3.943104e-01
res.aov <- anova_test(
  data = py$pdl_cr, dv = IC, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F        p p<.05   ges
## 1             chemo_regimen   1  72  7.410 8.00e-03     * 0.056
## 2               sample_type   1  72 22.106 1.21e-05     * 0.116
## 3 chemo_regimen:sample_type   1  72  1.398 2.41e-01       0.008
cys_den = pdl1_df3_log[pdl1_df3_log.sample_type == 'Cystectomy'].set_index('pt_number')
tur_den = pdl1_df3_log[pdl1_df3_log.sample_type == 'TUR'].set_index('pt_number')
common_in = cys_den.index.intersection(tur_den.index)
delta_den = cys_den.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1) - tur_den.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1)
delta_den = delta_den.assign(pt_number = delta_den.index)
delta_den = delta_den.assign(chemo_regimen = delta_den.pt_number.map(pt_reg)).assign(platinum_type = delta_den.pt_number.map(platinum_type_dict))

distances

log_dist3_cr = all_med_log[all_med_log.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]
log_dist3_ptt = all_med_log[all_med_log.platinum_type.isin(['Cisplatin', 'Carboplatin'])]
py$log_dist3_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(CD8_Tcell_to_tumor_median)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable                  statistic         p
##   <chr>       <chr>                <chr>                         <dbl>     <dbl>
## 1 Cystectomy  MVAC                 CD8_Tcell_to_tumor_median     0.946 0.120    
## 2 Cystectomy  gemcitabine_platinum CD8_Tcell_to_tumor_median     0.963 0.122    
## 3 TUR         MVAC                 CD8_Tcell_to_tumor_median     0.779 0.0000217
## 4 TUR         gemcitabine_platinum CD8_Tcell_to_tumor_median     0.951 0.0417
py$log_dist3_ptt %>%
  group_by(sample_type, platinum_type) %>%
  shapiro_test(CD8_Tcell_to_tumor_median)
## # A tibble: 4 × 5
##   sample_type platinum_type variable                  statistic           p
##   <chr>       <chr>         <chr>                         <dbl>       <dbl>
## 1 Cystectomy  Carboplatin   CD8_Tcell_to_tumor_median     0.936 0.574      
## 2 Cystectomy  Cisplatin     CD8_Tcell_to_tumor_median     0.956 0.0168     
## 3 TUR         Carboplatin   CD8_Tcell_to_tumor_median     0.912 0.366      
## 4 TUR         Cisplatin     CD8_Tcell_to_tumor_median     0.839 0.000000407
py$log_dist3_cr %>%
  group_by(sample_type) %>%
  levene_test(CD8_Tcell_to_tumor_median ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    78    0.0740 0.786
## 2 TUR             1    78    1.78   0.187
box_m(py$log_dist3_cr[, "CD8_Tcell_to_tumor_median", drop = FALSE], py$log_dist3_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      3.85  0.0497         1 Box's M-test for Homogeneity of Covariance Matric…
f1.ld.f1(y = py$log_dist3_cr$CD8_Tcell_to_tumor_median, time = py$log_dist3_cr$sample_type, group = py$log_dist3_cr$chemo_regimen, subject = py$log_dist3_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##              Statistic df      p-value
## Group       0.04884087  1 8.250927e-01
## Time       26.11852161  1 3.210892e-07
## Group:Time  0.04143483  1 8.387009e-01
f1.ld.f1(y = py$log_dist3_ptt$CD8_Tcell_to_tumor_median, time = py$log_dist3_ptt$sample_type, group = py$log_dist3_ptt$platinum_type, subject = py$log_dist3_ptt$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('Cisplatin', 'Carboplatin'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   Cisplatin Carboplatin 
##  If the order is not correct, specify the correct order in time.order or group.order.

##              Statistic df      p-value
## Group       0.04047519  1 8.405543e-01
## Time       22.55412414  1 2.043055e-06
## Group:Time  0.45706426  1 4.989992e-01
res.aov <- anova_test(
  data = py$log_dist3_cr, dv = CD8_Tcell_to_tumor_median, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F        p p<.05      ges
## 1             chemo_regimen   1  78  0.035 8.52e-01       0.000323
## 2               sample_type   1  78 17.819 6.51e-05     * 0.061000
## 3 chemo_regimen:sample_type   1  78  0.212 6.47e-01       0.000769
py$log_dist3_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(Macrophage_to_tumor_median)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable                   statistic        p
##   <chr>       <chr>                <chr>                          <dbl>    <dbl>
## 1 Cystectomy  MVAC                 Macrophage_to_tumor_median     0.846  4.10e-4
## 2 Cystectomy  gemcitabine_platinum Macrophage_to_tumor_median     0.931  6.47e-3
## 3 TUR         MVAC                 Macrophage_to_tumor_median     0.804  6.12e-5
## 4 TUR         gemcitabine_platinum Macrophage_to_tumor_median     0.855  2.54e-5
py$log_dist3_ptt %>%
  group_by(sample_type, platinum_type) %>%
  shapiro_test(Macrophage_to_tumor_median)
## # A tibble: 4 × 5
##   sample_type platinum_type variable                   statistic            p
##   <chr>       <chr>         <chr>                          <dbl>        <dbl>
## 1 Cystectomy  Carboplatin   Macrophage_to_tumor_median     0.959 0.796       
## 2 Cystectomy  Cisplatin     Macrophage_to_tumor_median     0.908 0.0000983   
## 3 TUR         Carboplatin   Macrophage_to_tumor_median     0.789 0.0220      
## 4 TUR         Cisplatin     Macrophage_to_tumor_median     0.804 0.0000000437
py$log_dist3_cr %>%
  group_by(sample_type) %>%
  levene_test(Macrophage_to_tumor_median ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    78     0.110 0.741
## 2 TUR             1    78     2.13  0.149
box_m(py$log_dist3_cr[, "Macrophage_to_tumor_median", drop = FALSE], py$log_dist3_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      4.33  0.0375         1 Box's M-test for Homogeneity of Covariance Matric…
f1.ld.f1(y = py$log_dist3_cr$Macrophage_to_tumor_median, time = py$log_dist3_cr$sample_type, group = py$log_dist3_cr$chemo_regimen, subject = py$log_dist3_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##             Statistic df      p-value
## Group       0.8742298  1 3.497870e-01
## Time       23.2285838  1 1.438434e-06
## Group:Time  0.1916301  1 6.615633e-01
f1.ld.f1(y = py$log_dist3_ptt$Macrophage_to_tumor_median, time = py$log_dist3_ptt$sample_type, group = py$log_dist3_ptt$platinum_type, subject = py$log_dist3_ptt$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('Cisplatin', 'Carboplatin'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   Cisplatin Carboplatin 
##  If the order is not correct, specify the correct order in time.order or group.order.

##              Statistic df      p-value
## Group       0.07567483  1 7.832466e-01
## Time       35.48936404  1 2.564470e-09
## Group:Time  2.54512061  1 1.106355e-01
res.aov <- anova_test(
  data = py$log_dist3_cr, dv = Macrophage_to_tumor_median, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F        p p<.05      ges
## 1             chemo_regimen   1  78  1.185 2.80e-01       0.010000
## 2               sample_type   1  78 18.644 4.58e-05     * 0.080000
## 3 chemo_regimen:sample_type   1  78  0.177 6.75e-01       0.000828